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^ ■ Abstract 

^ . An efficient integral equation based solver is constructed for the electrostatic problem on domains 

with cuboidal inclusions. It can be used to compute the polarizability of a dielectric cube in 
a dielectric background medium at virtually every permittivity ratio for which it exists. For 
example, polarizabilities accurate to between five and ten digits are obtained (as complex limits) 
for negative permittivity ratios in minutes on a standard workstation. In passing, the capacitance 
of the unit cube is determined with unprecedented accuracy. With full rigor, we develop a natural 
mathematical framework suited for the study of the polarizability of Lipschitz domains. Several 
aspects of polarizabilities and their representing measures are clarified, including limiting behavior 
both when approaching the support of the measure and when deforming smooth domains into a 

Q ■ non-smooth domain. The success of the mathematical theory is achieved through symmetrization 

c/3 , arguments for layer potentials. 
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Qs . The determination of polarizabilities and capacitances of inclusions of various shapes has a long 

0^ \ history in computational electromagnetics. Inclusions with smooth surfaces are, by now, rather 

standard to treat. When surfaces are non-smooth, however, the situation is different. Numerical 
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1. Introduction 



^J^ ' solvers can run into problems related to stability and resolution. Particularly so in three dimensions 

*vj . and for certain permittivity combinations. Solutions may not converge or results could be hard to 



interpret. See [4J, |45| . |51[ and references therein. The situation on the theoretical side is similar. 
When, and in what sense do solutions exist? Such questions are in the mainstream of contemporary 
research in harmonic analysis. Coincidentally, also in applied physics (plasmonics) there is a 
growing interest in solving electrostatic problems on domains with structural discontinuities and a 



X 

_Cp_' concern about the sufficiency of available solvers J19l. 154]. 



This paper addresses several fundamental issues related to the problems just mentioned. We 
construct a stable solver for the polarizability and capacitance of a cube based on an integral 
equation using the adjoint of the double layer potential. We compute solutions of unprecedented 
accuracy and interpret the results within a rigorous mathematical framework. The reason for 
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working with a cube are twofold. First, the cube has the advantage that its geometric difficulties 
are concentrated to edges and corners, since its faces are flat. Integral equation techniques, which 
often excel for boundary value problems in two dimensions, typically suffer from loss of accuracy in 
the discretization of weakly singular integral operators on curved surfaces in three dimensions. Here 
we need not worry about that. Secondly, cubes are actually common in plasmonic applications. 

In the purely theoretical sections we begin by collecting a number of results and recent advances 
in the theory of layer potentials associated with the Laplacian in Lipschitz domains. The most 
obvious reason for this is that the invertibility study of layer potentials leads to the solution of 
the boundary value problem implicit in the definition of polarizability, and is as such the basis 
for both the mathematical and numerical aspects of this paper. Furthermore, the properties of 
the polarizability for a non-smooth domain such as a cube are quite subtle, and it is our ambition 
to provide a solid theoretical foundation for the problem at hand, giving a careful and detailed 
exposition of a mathematical framework that clarifies a number of points. 

Since the double layer potential is not self-adjoint in the L^-pairing, we develop certain sym- 
metrization techniques for it, in particular extending the work of Khavinson, Putinar and Shapiro 



311 ] to the case of a non-smooth domain. These techniques are used to prove the unique existence of 
the polarizability itself for a Lipschitz domain, as well as of a corresponding representing measure 
[ig ]. We present a thorough discussion of the smooth case, the limiting behavior in passing from 
the smooth to the non-smooth case, and ultimately the general case. Concerning the last point, 
a condition ensuring that the representing measure has no singular part is given, and it is proven 
that in the support of the absolutely continuous part of the measure, the polarizability can not be 
given a direct interpretation in terms of a potential with finite energy solving the related boundary 
value problem. 

The paper is organized as follows: Section [2] formulates the electrostatic problem and defines the 
polarizability. Existence issues and representations are reviewed in Section [3j For ease of reading, 
rigorous statements and proofs are deferred to Sections H] and [5j The capacitance is discussed in 
Section [6l Section [7] reviews the state of the art with regard to numerical schemes. Section [8] gives 
a necessary background to the present solver. New development takes place in Section [9j The last 
sections contain numerical examples performed in Matlab. Section [10] illustrates the effects of 
rounding corners and Section [11] is about the cube. 

The main conclusion of the paper is that, from a numerical viewpoint, it is an advantage to let 
cubes have sharp edges and corners as opposed to the common practice of rounding them slightly. 
Furthermore, the representing measure for the polarizability of the cube is determined, and a new 
benchmark for the capacitance of the unit cube is established. 

2. The electrostatic problem and the polarizability 

Let a domain V, an inclusion with surface S and permittivity 62, be embedded in an infinite 
space. The exterior to the closure of V is denoted E and has permittivity ei. Let z/^ be the exterior 
unit normal of S at position r. 

We seek a potential U{r), continuous in EL) S L)V, which satisfies the electrostatic equation 

AC/(r) = 0, r£EUV, (1) 

subject to the boundary conditions on the limits of normal derivatives 

^^^U'^\r) = e2^U'-\r) (2) 



and behavior at infinity 

lim VU{r) = e. (3) 

1 — >oo 

Here superscripts ext and int denote limits from the exterior or interior of 5, respectively, and e is 
an applied unit field. Eqs. ([I]), ©i dSl) constitute a partial differential equation formulation of the 
electrostatic problem. Proposition 15.11 gives a strict interpretation of what it means for a potential 
U{r) to solve this problem, in particular expressing ([2]) in a distribution sense. 

For the construction of solutions to ([T]) , ([2]) , ([3]) we make use of fundamental solutions to the 
Laplace equation in two and three dimensions 

G(r,r') = loglr — r'l and G(r,r') = -—- , (4) 

27r 47r \r — r'\ 

and represent U{r) in terms of a single layer density p{r) as 

U{r) = e • r + / G{r, r')p{r') dar' , (5) 

Js 

where da is an element of surface area. 

The representation ([5]) satisfies ([1]) and ([3]). Its insertion in ([2]) gives the integral equation for 
p{r) 

pir) + 2X [ ^G{r,r')p{r')dar' = -2X{e-Ur) , reS, (6) 

where the parameter 

X='-^^^. (7) 

£2 + 61 ^ ^ 

The polarizability tensor of V can be defined in terms of an integral over a polarization field. 
When V features sufficient symmetry, such as the octahedral symmetry of the cube, the polariz- 
ability is isotropic and reduces to a scalar a(ei,e2), see [48|, which can be determined via 

(8) 



(e2-ei) / VU{r)dr = a(ei,e2)e, 
Jv 



where dr is a volume element. Using integration by parts in ([8]) it is possible to express a(ei,e2) 
as an integral over p[r) 



a(ei, £2) = -ei / p{r) (e • r) dar ■ (9) 

Js 

More generally, the components of the polarizability tensor can be recovered via integrals similar 
to that in ([9]) using different applied fields e, but in what follows we tacitly assume that V is 
sufficiently symmetric as to allow for a scalar valued 0(61,62). 

3. Theory — overview 

For what surface shapes S and permittivities 61, 62 does the electrostatic problem have a 



solution? Starting with ([6]) and partly following 37|, this section sketches the derivation of some 



important existence results. It also motivates an integral representation formula for 0(61,62) and 
two sum rules which are used for validation in our numerical experiments. 



3.1. Existence of solutions for smooth S 

Let us rewrite (161) in the abbreviated form 



{I + XK) pir) = Xgir) , (10) 



where / is the identity. If S is smooth, then (|10p is a Fredholm second kind integral equation 
with a compact, non-self adjoint, integral operator K whose spectrum is discrete and accumulates 
at zero. Let K and its adjoint K* (the double layer potential) have eigenvectors (pi and ipi with 
corresponding eigenvalues Zj. All eigenvalues are real and bounded by one in modulus. Non-zero 
eigenvalues have finite multiplicities. Normalizing 



%lji{r)(t>j{r)dar = Sij, (11) 

s 

the kernel of K can be written 

K{r,r') = Y,ZiMr)'^^. (12) 



See, further, the discussion in Section 1 

Let us introduce a new variable z and a scaled polarizability a{z) as 

z = -l/X, (13) 

a(^) = %f^, (14) 

where \V\ is the volume of V. Then (114p . with Q, can be written in the abbreviated form 



a{z) = / h{r)p{r)dar. (15) 

Js 

The relation (fT3|) allows us to use the parameter A or its negative reciprocal z, depending on what 
is most convenient in a given situation. 
In terms of the quantities 



h{r)(j)i{r) dar and Vi = %lji{r)g{r)dar , (16) 

5 Js 



and using (|12p to construct the resolvent of (jlUp . one can write 

i 

and 



«(-) = E7^' (18) 



see Theorem 15.61 This suggests that neither U{r) nor a(z) exists for z = Zi when UiVi 7^ 0. There 
is an electrostatic resonance or plasmon at Zj. 

For ease of interpretation, the sum in (|18p can be considered as taken over distinct eigenvalues 
and with UiVi, for a degenerate eigenvalue, being the sum of all residues belonging to that eigenvalue. 

4 



Then all UiVi are non-negative and plasmons can be classified as bright or dark depending on whether 
UiVi > or not [5J]. When 5" is a circle, there are only two eigenvalues: zi = —1 which is simple 
and corresponds to a dark plasnion and Z2 = which has infinite multiplicity and corresponds to 
a bright plasnion. When 5 is a sphere, the eigenvalues are Zi = 1/(1 — 2i). The multiplicity of Zj 
is 2i — 1. The only bright plasmon is associated with Z2- 

For later reference we observe that the sum of all residues is 

^UiVi= f h{r)g{r)dar = 2. (19) 



When the polarizability is isotropic one can, using techniques from [151], also derive a weighted sum 
rule 

Y^ ZiUiVi = j I h{r)K{r, r')g{r') da', da, = 2{2/d - 1) , (20) 

where (i = 2, 3 is the dimension. 

3.2. Existence of solutions for non-smooth S 

If S is gradually transformed from a smooth surface into a non-smooth surface, eigenvalues 
Zi travel and occupy a certain subset of the interval [—1,1] ever more densely. When S ceases 
to be smooth, K is no longer compact with discrete eigenvalues. Rather, K has a continuous 
spectrum which on a certain function space coincides with the aforementioned subset, accompanied 
by discrete values. Disregarding the discrete spectrum, which for squares and cubes turns out to 
correspond to dark plasmons, the sum (|18p assumes a limit 



a(z)= / ^^= / ^^^Mi^, (21) 



where the measure ^{x) is real and non-negative and a^ = {x : fJ-'{x) > 0}. Here we have ignored 
the possible presence of a singular spectrum. For further details and a condition that serves to 
exclude this complication, see Theorems 15.21 and Section 15.21 The sum rules ()19p and ()20p assume 
the forms 

/ ^'(3;)dx = 2, (22) 

/ V(x) dx = 2(2/d - 1) . (23) 

Jr 

The numerical results in Sections [10] and [11] suggest that both the square and the cube have 
a^ equal to a single, possibly punctured, interval (a, 6) C [—1,1]. The square has a = —0.5 and 
b = 0.5, consistent with the exact computations of the spectral radius for a square found in [35| 



and [53|. The cube has a « —0.694526 and b = 0.5. For later reference we let u^sq denote a^ of 
the square and u^cu denote o"^ of the cube. 

For a large class of non-smooth S, the potential U{r) exists when z stays away from a certain 
compact set L : (T^j C L C [— 1, 1]. Furthermore, a{z) has a limit, 

a'^{x)= lim a{x + iy), (24) 



as z = X + iy approaches x from the upper half-plane for almost all x G M. See Theorem 
and Section 15.21 It is important in this context and when x E o"^ not to interpret a'^{x) as a 
polarizability corresponding to a meaningful solution U{r) for a negative permittivity ratio £2/61 = 
{x — l)/{x + 1). On the contrary, Theorem 15.81 states that there is no U{r) with finite energy 
solving (P), (121), dSI) when z = x £ a^. Therefore, any attempt to solve the electrostatic problem 
directly at a point z £ a^ is bound to fail. 

3.3. The limit polarizability a^(x) and its relation to fi'{x) 

This paper aims at constructing an efficient scheme for computing a{z) of a cube at all z for 
which this quantity exists. Still, in our numerical experiments we only compute the limit a~^{x) 
of (|24p for X € [—1, 1]. The reason for this is that the computation of a{z) is hardest for z close to 
a^ C [—1,1]. Accurate results for a+(x) therefore indicate a robust scheme. Furthermore, there is 
a simple connection between a^{x) and n'{x). Using jump relations for Cauchy-type integrals one 
can show from ()2ip that 

^'(x) = 9{a+(x)}/7r, xeR. (25) 

Knowledge of n'{x) for non-smooth S is of great interest in theoretical materials science. Closed 
form expressions seem to be out of reach, however, except for a famous example in a periodic two- 
dimensional setting [lO|, |40|. As for numerics, merely determining a^ is a challenge 4J, |51|. To 



the authors' knowledge, a^ is not known for any S exhibiting corners in three dimensions. The 



accurate determination of fi'{x) is even harder 30(]. Studying how psp evolves as a smooth S 



becomes non-smooth is not an efficient method. Eigenvalue problems are costly to solve. The 
discretization of K on surface portions of high curvature is problematic. Conditioning is also an 
issue and details of the mapping UiVi — )• /u'(x) need to be worked out. It is desirable to find /i'(x) 
in a more direct way and ()25p offers precisely this. Obtaining fi'{x) is a subproblem of computing 

a+(x) for X £ [—1, 1]. 

4. Theory — preliminaries 

Let V C M*^, d > 2, be an open and bounded set that is Lipschitz, in the sense that its boundary 
S = dV is connected and locally the graph of a Lipschitz function in some basis. For a more precise 



definition of this concept, see for example 50|. To avoid a certain technicality we will also assume 
that V is star-like, meaning that there exists an r^ £ V such that the line segments between tq 
and every other point r £ V are contained in V. In the applications of this paper, V will take on 
the role of the square in R^, the cube in M^, or a set with smooth boundary approximating either 
of the two. As before we will denote E = V . 

In this section we first record a number of results about the single and double layer potentials 
associated with the Laplacian on V, to then introduce the mathematical framework in which we 
will study the boundary value problem given by ([T]), ([2]), ([3]). Actually, in this section and the next, 
we will develop the theory only for d > 3. The two-dimensional case contains several anomalies 
in relation to the higher-dimensional theory, and it is for the sake of clarity and brevity that we 
exclude it. We shall indicate some of the differences as we progress, but we note here that the 
main results about the polarizability a{z) remain true also for d = 2. 

While L'^{S) is a natural domain for the operator K, we will primarily focus on the action of 
K on certain Sobolev spaces H'^. There is good reason for this. For one, K is not self-adjoint as 
an operator on L'^{S), or even normal, so that the spectral theorem can not be directly applied. 
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We will therefore develop certain symmetrization techniques, and these demand that we consider 
K on fractional Sobolev spaces. A second, related reason, is that the -L^-spectrum of K is no 
longer contained in the real line when S fails to be smooth, see I. Mitrea [,42|. We shall see that 
considering K ona Sobolev space amends this problem. We also note here that by the X-spectrum 
of a bounded operator T on a Hilbert space X, T : X — )• X, we always mean the set 

Spec(T, X) = {z G C : K — z is not bijective on X}. 

When z ^ Spec(T, X) it is a consequence of the closed graph theorem that K — z has a bounded 
inverse {K — z)~^ : X ^ X. 

For s = we have simply that H^{V) = L'^{V) and H^{S) = L'^{S). For s = 1, H^{V) is the 
Hilbert space of distributions u such that u and dx u, 1 < j < n, are members of L'^{V). The norm 
is given by 

ll^ll//l(y) = ll^llL2(y) + ||Vn||^2(y). 

H^{S) can be defined similarly using the almost everywhere defined tangential vectors of S, see 
for example Geymonat [iTI]. For < s < 1, H^ can be defined by real interpolation methods, but 
in our situation it can alternatively be characterized by a Besov type norm. That is, u G H'^{V) if 

II ii2 II ii2 f Wi^) - uir')]"^ , ^ f 

JVxV I' ' I 

u G H^(y) is then inductively defined for s = Si + Sf with Sj > 1 an integer and < sj < 1 by 
requiring that u G H^^ and d^u G H'^f for /3 = (/3i, . . . , /3rf) with Ylil^k = Si. Returning to the case 
< s < 1, we have that u G H''{S) if 

II ii2 II ii2 /" l^('') -■u(r')p , 

\mH'{S) = m\v^{S) + / \^_^,\d-i+2s d^- da,. < oo, 

JSxS I' ' I 

where a denotes Hausdorff measure on S. See Adams [V\ and Grisvard 20] for further information 
and the equivalence of various definitions in the Lipschitz setting. For s > we define H~^ as the 
dual space of i?* in the L^-pairing. More precisely, a distribution u lies in H~'^ if and only if 

ll^ll//-" = sup \{u,v)i2\ < oo. 

We shall also make use of Sobolev traces, which give us a way to assign boundary values to distri- 
butions in V. We will only require the classical Gagliardo result [la] which says that there uniquely 
exists a continuous, surjective linear operator Tr : H^{V) — t- H^'^{S) with right continuous inverse 
such that Tru = u\s for any u G C°°{V). There is also a corresponding trace from the exterior 
domain with the same properties, Trg : H^{E) — )• H^''^{S). 

The L^(S')-adjoint K* of the operator K is known as the double layer potential, given by the 
formula 

{K*u){r) = 2 [ -^G{r, r')u{r') da,/, u Ci L^{S), r £ S. (26) 

Js OVr' 

For d = 2,3 the Newtonian kernel G has already been defined in dH), and for d > 3 it is given by 



G{r,r') = ujd\r — r' 



l\2-d 



with a normalization constant oj^ chosen so that ArG{r, 0) = —6 in the sense of distributions. 
When 5 is a C'^-surface the kernel of K* is only weakly singular (for d = 2 there is no singularity 
at all present), and it is a standard matter to see that ()26p defines K* as a compact operator on 
H^{S) for < s < 1. The compactness of K* makes its spectral analysis considerably easier, and 
in this case it is well known that 

Spec(i^*,L2(S))c[-l,l), (27) 

see for example the results of Escauriaza, Fabes and Verchota [ll|] together with the fact that the 
spectrum of K* is real in the C'^-case. This latter point will be discussed further later on. 

Unfortunately, when S is only a Lipschitz surface, K* is no longer compact in general. In fact, 
when S" is a curvilinear polygon in two dimensions, I. Mitrea [4^ has shown that the L^-spectrum 
of K* consists of the union of certain solid "figure eights" in the complex plane, one for each 
(non-smooth) vertex of S, in addition to a finite number of real eigenvalues. In particular this 
applies when S* is a square in two dimensions, with only one figure eight present, since all angles 
are equal. The general situation is not as well understood, but when y C M'^ is convex, as it is in 
our situation, it is known that the spectral radius of K* on L'^(S) is 1, see Fabes, Sand and Seo 

0. 

To even define K* in the general Lipschitz setting, the integral in ()26p must be understood 
in an almost everywhere principal value sense. We remark, however, that when 5" is a curvilinear 
polyhedron, K*u{r) can be evaluated in the usual integral sense, except possibly when r belongs to 
an edge of S, and so it is not necessary to consider principal values in the main applications of this 
paper. Proving the boundedness of K* on L'^{S) was an accomplishment of Coifman, Mcintosh 
and Meyer [9[ in their study of singular integrals. The boundedness of K* as an operator on H^{S), 



< s < 1, also essentially follows from [9|, see for example Meyer [38|]. By duality we immediately 
obtain that K is bounded on H~'^{S), < s < 1. 

For u G L'^iS) one may of course also evaluate the integral (|26p in V U E to obtain a harmonic 
function. We denote 

{Du){r) = 2 / 7^— G(r, r')u{r') dd^, reViJE. 

J S '^Vr' 

Fabes, Mendez and M. Mitrea 112] prove that for < s < 1, D : i?*(5) -> H''^'^/'^{V) is bounded. 
The single layer potential of u, defined in all of W^, is given by 



{Su){r) = 2 / G{r,r')u{r')dar', u G L\S), r G M*^. 
Js 

The kernel G is only weakly singular when 5" is a Lipschitz surface, so that there is no issue in 



defining this integral operator. In fact, S has smoothening properties. D. Mitrea 4l|] shows that 
for < s < 1, 5 : H~^{S) — )• H^~^{S) is a bicontinuous isomorphism and in [1^ it is proven that 

3 

S is bounded as a map S : H^^{S) — )• H2~'^(y) for < s < 1. It is clear that S is self-adjoint in 
the L^(5)-pairing and that Su is harmonic inVUE. 

At this stage, a peculiarity of the case d = 2 appears. In any dimension, there exists uniquely a 
function uq £ L'^{S) such that {I + K)uq = and JgUQ da = 1, and one can show that 5no|y = c is 
constant. In higher dimensions this constant can never be zero, but for d = 2 there exist domains 
such that c = 0. When this occurs S clearly fails to be injective, and its range is also affected. On 
the other hand, if c = for a particular domain V, any non-trivial dilation of V will give a domain 



with c 7^ and the properties in the previous paragraph may be proven to hold for the dilated 
domain, at least for s = 1/2, which will tm'n out to be the important case for us. See Verchota 
[501 1 for details. Note that since our object of interest, the polar izability a{z), is scaling invariant, 
this anomaly of S for d = 2 presents no real obstacle. 

While the kernel of S is sufficiently nonsingular to immediately define a continuous function 
Su everywhere on R if u is for example bounded, similar statements are never true for the kernel 
of K* . In fact, the following jump formulas hold for a function u S L'^{S). 

d^S^^'^u = -u + Ku D'^'^u = -u + K*u (28) 

D^^t^ = u + K*u, 

where a superscript int or ext denotes taking a limit from the interior or exterior of S, respectively. 
In general the formulas are true in the sense of non-tangential convergence almost everywhere on 
S, see [50]. These jump relations explain why the boundary condition 1^ leads to the integral 
equation ([5]). In a moment we shall make a more precise statement about this. Before doing so, 
we need to show that the jump formulas hold in a certain trace sense. 

For this purpose, we will also need to consider the Hilbert space T-L{V)/C of harmonic functions 
V on V, modulo constants, with finite energy, 

ll^llw(y)/c = / iVul^dr <oo. 

Since this semi-norm annihilates constants we consider v and w-|-C, c G C, to be the same element. 
Note that 'H{V)/C is continuously contained in H^{V)/C by the classical Poincare inequality for V. 
Since V is assumed star-like, it is straightforward to use dilations in order to prove that functions 
which are harmonic and smooth in V are dense in H^{V)/C In fact, we introduced the hypothesis 
that V is star-like only to facilitate such density statements. 
The Dirichlet problem 

ve ?{{¥)/€ Tiv = u, 

is well-posed for initial data u G H^''^{S), see for example [12i]. Equivalently, Tr : 'H(y)/C -^ 
H^''^{S)/C is a bicontinuous isomorphism. Often we will simply denote Trv = v\s when it is clear 
what is meant. 

It is established in a paper by Hofmann, Mitrea and Taylor [28] that Green's formula 



{V(t>,Vilj)dr+ / (/)AV^dr= / 05^^ do" (29) 

V Jv Js 

continues to hold true for S Lipschitz and (j),ip G C°°{V). Since the functions in T-L{V)/C are 
harmonic, this shows that its scalar product satisfies 

(■u,w)^(y)/C = / {'Vv,Vw)dr = / vd„wda= / {duv)wda. 
Jv Js Js 



Initially these identities are valid only for smooth v and w, but as in [3l|] one can argue by duality 
and density to interpret the normal derivatives d/yV and d/^w as elements of H~^i'^{S) so that the 
equalities remain true. If we denote by Hq (5) the closed subspace of those u G H^^''^{S) 

9 



such that jguda = 0, the imphed duahty argument gives rise to a bicontinuous bijective operator 

d,^ : H^''^{S)/C — )• Hq (S) which should be understood as the normal derivative of the trace. 

It is important for our purposes to now repeat this construction for the exterior space %{E) 
of harmonic functions v va. E with finite energy norm and \m\r-^oov{r) = 0. We state this as a 
proposition. 

Proposition 4.1. The exterior trace is a bicontinuous isomorphism when considered as an operator 
TtE '■ T^{E) — )• H^''^{S). There is a corresponding bounded bijective operator d^ : H^''^{S) — )• 
H~^'^{S) satisfying 

{v,w)'^(^E) = / (V^;,Vu))dr = — / vdj^wda = — {d^v)w(ia. 
Je Js Js 

Proof. The statements about Tr^; follow by the well-posedness of the Dirichlet problem, see [12l |. 
The construction of d^ again follows along the lines of [3l|] . D 

Remark 4.1. When d = 2 the additional condition Jgudcr = is required to solve the exterior 
Dirichlet problem v G 7i{E), Tiev = u. This is also reflected in the kernel G{r,r') of the single 
layer potential. Note that G{r,r') ~ — ^log|r| as r — )• oo for d = 2, but linir-^oo G{r,r') = for 
d>2. 

We end this section with an interpretation of the jump relations ()28p within the just established 
framework. 

Proposition 4.2. Let u £ H^^^{S), then Du G %{¥)/€, Du G n{E) and 

TtDu = -u + K*u TteDu = u + K*u. 

Furthermore, let v G H-^/'^{S). Then Sv G n{V)/C, Sv G TiiE) and 

TrSv = TrESv = Sv\s d^Sv = v + Kv d^Sv = -v + Kv. 

Proof. Suppose first that u and v are smooth and harmonic on V . Then by applying Green's 
formula we find that 

Du[r) = S{duu){r) - 2u{r), r £V. (30) 

In particular, this shows that Du\y extends continuously to V ., and hence the jump relation Z)'"*u = 
—u + K*u must hold in trace sense. That is, Tr Du = —u + K*u. Since both sides of this equation 
are continuous maps of H^''^{S) by previously quoted results, it must hold for every u G H^''^{S). 
By the same reasoning we obtain TiSv = Sv\s for every v G H^^''^{S). 

To show that d^Sv = v + Kv we note that S{dyu) = u + K*u by PUJ) and the jump formula. 
The sought formula is the dual statement of this. More precisely, for any smooth harmonic ip we 
have 

{duSv,ip)L2(s) = {Sv,dyij)L2^s) = {v,S{d,yip))L2(s) 

= {v,^ + K*'ip)L2^s) = {v + Kv,ij)L2(s), 

which verifies that d^Sv = v + Kv for all v G H^^''^{S) by continuity and density. 

The exterior statements are dealt with similarly. D 
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5. Theory — results 

5.1. Existence of the measure /x 

We are now in a position to develop the symmetrization techniques that have been alluded to 
previously. Once these are in place, we can use the spectral theory of self-adjoint operators to 
prove that the (scaled) polarizability a{z) = \y\'^'^ , z = T^f^ G C, has a representing measure fi. 

We begin, however, by describing the sense in which the potential U will solve the boundary 
value problem given by ([1]), ([2]) and ([3]). Note that the following proposition furthermore expresses 
the fact that if looking for a potential such that U{r) — e ■ r has finite energy, then H^^''^{S) is 
exactly the right space to find the corresponding density distribution p. 

Proposition 5.1. Let p € H^^''^{S) be such that 1^ holds, i.e. {K — z)p = g, where g{r) = 
— 2(e ■ Ur)- Let 

U{r)=e-r + -Sp{r), r G M^. (31) 

Then U G 'H{V)/C, U - e ■ r e TiiE), TrU = TteU, lim^^oo Vf/ = e and U satisfies ^ in the 
sense that 

ei (af ([/ - e • r) + d^{e ■ r)) = es^^f/. (32) 

The converse is also true. That is, if U satisfies the above properties, then there exists a p G 
H-^/^{S) such that (^ and © hold. 

Proof. This is a consequence of Proposition 14.21 the well-posedness of the interior and exterior 
Dirichlet problems and the bijectivity of S : H^^/^{S) -)> H^/^{S). D 

— 1/2 

Remark 5.1. When z ^ —1, the hypothesis that {K — z)p = g implies that p G Hq {S). 

— 1/2 

This seen by taking into account that g = —2d,y{e ■ r) and d^Sp both belong to Hq [S) in the 
computation 

z pda = Kp — gda = / Kp da = dySp — pda = — / p da. 

Js Js Js Js Js 

This is of importance for the case d = 2 (cf. Remark 14. ip . 

In the sequel we shall denote p = pz and U = Uz to indicate their dependence on z. Under 
the hypothesis of the preceding proposition, we can, due to the assumption of isotropy, express the 
scaled polarizability as 

"(^) = Tf^ I '^Uz{r)-edr = %-^ f {d,Uz){r){e ■ r) da,. 
|i/|€i Jv r Ki Js 

e2-ei /", ,1 



iT.| {e-Ur + -{pz + Kpz){r)){e ■ r) da,. 

€2- ei z + 1 



T^i o , Pz{r){e ■ r) dar = / Pzhda, 

>^|ei ^ Js Js 

where h{r) = —{e-r)/\V\. Since {K — z)pz = g, this shows that the analysis of a is closely related 
to the spectral theory of K. 

The scalar product on T-L{V)/C is one of the keys to understanding the spectral theory of K 
and K*, since K* is self-adjoint in the ^(y)/C-pairing. To be precise, the above shows that any 
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element v G H^''^{S)/C can also be considered as an element v G 'H{V)/C in a bicontinuous way. 
We are therefore justified in letting H^''^{S)/C inherit its scalar product from 'H{V)/C, 

{v,w)hi/2(^s)/c = {v,w)n(v)/c- 

Since K* maps constants onto constants (cf. ()30p ). we may consider K* as a bounded map on 
H^/^{S)/C. Let v,we H^/^{S)/C. By the fact that K*v = S{dyv) - v ii then holds that 



{K*v,w)hi/2(^S)/'C = l^ii'S{duv) - v){d^w)da 

vdu{S{dyw) -w)da = {v,K*w) ^i/2^syQ. 



s 



Theorem 5.2. There exists a compact set L C M and a positive Borel measure fi with total mass 
2 and compact support contained in L, such that for z ^ C, z ^ L, [K — z)p = g has a unique 



— 1/2 

solution pz G Hq (S) and 



a{z) = fM^. (33) 

Jr X- z 

fj, is unique in the class of compactly supported finite Borel measures such that (j33p holds for all 

z e C+ = {z : '=jz > 0}. 

Proof. Since K* is bounded and self-adjoint on H^''^{S)/C it has by the spectral theorem a 
corresponding projection- valued spectral measure £ with support in the spectrum of K*. Let 
L = Spec-ff*. £ is characterized by the fact that for every bounded Borel- measurable function / 
on L, it holds that 



f{K*) = j f{x)d£{x) 



(34) 



For any z ^ L note that K* — z \s invertible on H^''^{S)/C, and hence K — z is invertible on the 
dual space H~^^'^{S). Since h G H^/^{S) and g = 2\V\dyh G H^^^^{S) we have 



aiz) = {{K - z)-'g, h)L2(s) = {9, {K* - zy'h)L2(s) 

= 2\V\{d,h, {K* - z)-'h)L2^s) = 2\V\{h, {K* - z)-^M^i/2(5)/c 

^(S)/C- 



2\V\{{K*-z)-'h,h)j,^;,. 



With p{R) = 2\V\{£{R)h,h)jji/2/s)/c ^^^ ^^y Borel set i? C M we obtain by applying ^^ with 
f{x) = l/{x — z) the desired formula 

f dp{x) 



X — z 



The constructed measure p is positive since fJ-{R) = 2\V\{£{R)h,h) = 2\V\{£{R)h,£{R)h) = 
2\V\\\£{R)h\\'^ > 0, and since £(W) is the identity operator, its total mass is 

||/.|| = 2|y|||/i||2 = 2-^ /|epdr = 2. 
\^ \ Jv 

Equation (|33|) expresses that a is the Cauchy transform of p, a{z) = K,^{z). See Cima, Matheson 
and Ross [8| for an excellent survey of the Cauchy transform, written for measures with support 
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on the unit circle, but all results can be transformed to results about measures on R through 
standard conformal mapping techniques. See Koosis [33] for an explanation of this latter point, 
as well as results stated directly for the real line. By the classical F. and M. Riesz theorem, any 
other measure fi* supported on M and such that /C^* (z) = a{z) = JC^{z) for z £ C+ must be of the 
form d^* = d/i + udx, where u £ 7i^{C+) is given by the boundary values of an analytic Hardy 
space-function in the upper half-plane. Since such a function u can never be compactly supported 
unless it is identically zero, we obtain the uniqueness part of the theorem. D 

Remark 5.2. The fact that ||/i|| = 2 is sum rule ([22]). Note also that L = Spec{K*,H^/'^ /C) = 
Spec{K,HQ ) only differs from Spec{K* , H^'^) = Spec{K,H~^''^) by the point z = —1, which is 



in the latter spectrum but not in the former, see 4l|. 



The considerations leading up to the previous theorem are quite similar in spirit to those of 
Bergman [4]. Bergman considers a different potential operator which is symmetric under the inner 
product Jy{Vv, Vw) dr, although the arguments presented are somewhat incomplete since a space 
of functions belonging to this inner product is not identified. 

Before discussing the specific features of a and ^ we shall gain some further insight into the 
structure of K and K* by investigating a different symmetrization approach, expounded upon 
in the case when 5 is a C^-surface by Khavinson, Putinar and Shapiro [31]]. The starting point 
is Plemelj's symmetrization principle, which says that SK = K*S on L'^{S) and continues to 
hold true in our situation with the same proof as in [3]|. This operator equality amounts to the 
statement that K is self-adjoint under the inner product {Su,v)i2(^sy Note that 5 is a strictly 
positive operator on L'^{S), so that the form {u,v) — >• {Su,v)i2(^g-^ is strictly positive definite. 

Instead of working with the completion of L^ ( 5) under {Su , u) 2,2 (5) , we shall follow the approach 



of 31[ and introduce an operator-theoretic formalism to express the symmetrization of K. Recall 
that K : L'^{S) — )• L'^{S) is compact when 5 is C^, and so its spectrum consists of the point 
and a sequence (zj) of non-zero eigenvalues tending to zero, every eigenspace Hz{K) = ker(i^ — z) 
having finite dimension when z 7^ 0. 



Theorem 5.3 ([3]|). Suppose that S is a C'^-surface. Then there exists a self-adjoint compact op- 
erator A : L'^{S) -^ L2(5) such that AVS = VSK on L'^{S). When z / 0, \/S : H-,{K) ^ Hz{A) 
and VS : Hz{A) — )• Hz{K*) are isomorphisms of the indicated eigenspaces, and the eigenvectors of 
K* (including those for z = 0) span L?'{S). In particular, the LF'^S)- spectrum of K* is real. 

From the point of view of Khavinson et al., the proof of this theorem parallels the symmetriza- 
tion theory of Krein [34l | , w hich mainly concerns compact operators. However, the theory of Hassi, 
Sebestyen and de Snoo [2l|] assures us of the existence of A even when K is only bounded, that is, 
when S is only Lipschitz. The key fact is still that SK = K*S. 

Proposition 5.4. There exists a hounded self-adjoint operator A : L^{S) — )• L'^{S) such that 
AVS = VSK onL^{S). 

Proof. Given the existence of a bounded operator A on Lp'{S) satisfying A\/S = VSK, we deduce 
from the equation \'SA\'o = SK = K*S = \/SA*\/S, and the injectivity and dense range of S, 
that A must in fact be self-adjoint. D 



In fact, we have the following result, contained in 3l|, Proposition 1] with a proof that carries 
over to our setting. 
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Proposition 5.5. yS : L'^{S) — t- H^''^{S) is a bicontinuous isomorphism. Dually, yS also extends 
to a bicontinuous isomorphism \/S : H~^''^{S) — )• L^{S). 



«(^) = /r^- Since 



The existence of A gives us another way to derive the existence of a measure /i such that 

[K-z]-^ = VS~\A-z)-^VS (35) 

as an inverse on H^^''^{S), we can write 

a{z) = {{K - z)-'g, h)L2^s) = {{A - z)-'VSg, Vs'\)l2^s), (36) 

where, clearly, \/Sg,\/S h E L?{S). If £a is the spectral measure of A, we therefore obtain a 
measure with the desired property by letting /i(i?) = {£A{R)VSg, VS /i)l2(5) for i? C M. In this 
way we immediately recover all the conclusions of Theorem 15.21 except for the positivity of /i. 

To see the positivity, recall for u G H^''^{S) that S{dyu) = u + K*u and consider the compu- 
tation 

VS{d^u) = VS~^{u + K*u) = {I + A)VS^^u. 

In view of the identity g = 2\V\duh we find that 

a{z) = 2\V\{{A + I){A-z)-^^'^h,Vs'\)L2isy 

Clearly, it follows that /x>Oif^ + />0, that is, if 

Spec(K, H~^/'^{S)) = Spec(^, L^{S)) C [-1, oo). (37) 

In the C^-case we know from (j27p and the compactness of K that Spec(/f, iJ~^'^) C [—1, 1]. We 
shall see in Theorem 15.71 that this continues to be true when S is merely a Lipschitz surface. 

5.2. Properties of the polarizability a{z) and the measure fi 

In the final part of this section we will employ all the tools introduced thus far, in order to 
understand the behavior and specific features of the polarizability a{z) and its related measure 
fi. We begin by discussing the case of a smooth surface S, giving a rigorous treatment of many 



of the formulas and ideas that appear in the paper 37[ by Mayergoyz, Predkin and Zhang. From 
there we discuss the idea of approximating a non-smooth surface by a sequence of smooth surfaces, 
proving that the corresponding representing measures converge in a weak-star sense. We conclude 
with an analysis of the measure /U in the general non-smooth case, in particular giving a condition 
which guarantees that it does not have a singular part. Furthermore, we show that even though 
the polarizability a~^{x) exists in a limit sense almost everywhere x G M, the potential Ux does not 
exist at points where fJ,'{x) > 0, neither as a solution of the boundary value problem given by ([1]), 
(121) and (l3|), nor in a limit sense. 

Suppose now that 5" is a C^-surface, so that the operator A is compact. Let (f^) be an 
orthonormal basis for ker A, and let (fl) be eigenvectors of A corresponding to the non-zero 
eigenvalues (zj), repeated according to multiplicity, so that (/j) = (/j?) U (//) is an orthonormal 
basis for L'^{S). By Theorem 15.31 K,K* : L^{S) — )• L'^{S) have the same non-zero eigenvalues (zj), 

and corresponding (non-orthonormal) eigenvectors are obtained as (pj = v5 fl and ipj = \/Sfl, 
respectively. Concerning zero eigenvectors, ijj^ = ySff € H^''^{S) are clearly in the kernel of 
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K* , but 4>^ = v5 ff are in general elements of H^^''^{S) and only zero eigenvectors of K when 
considered as an operator on said space H~^''^{S). In particular, we note that the H^^''^{S)- 
eigenvectors of K span the whole space H~^''^{S). With this in mind, we shall denote ((/>j) = 
((/>?) U {(j)]) and (^j) = (^j°) U {i^}), the indexing arranged appropriately so that 

{(pi,i^j)LHS) =^ij- 

As a final notational detail, we shall let (zj) denote the full sequence of eigenvalues including zeros, 
so that K(j)i = Zi<j)i and K*^pi = Ziipi for every i. 

Based on the spectral decomposition of A we obtain for u £ Lp'{S) that 

^fSKu = A^fSu = ^ Zi{VSu, fi)L^{s)fi, 

i 

with convergence in L?'{S). Equivalently, we have for u G H^^''^{S) the expansion 

Ku = J2 ^ii'^^ i'i)L^{S)4>u (38) 

i 

with convergence in H~^''^{S). This is the formal interpretation of (|12|) . In this framework it is 
now easy to furthermore justify (fT7|) . (fTHjl and (fT9]) . We state this as a theorem. 



Theorem 5.6. Let S be a C'^ -surf ace. Then ()38p holds with K considered as an operator on 
H~^''^{S). Furthermore, let Ui = {(j)i,h)i^2/g\ and Vi = {g,'4)i) 1,21 gy Then 

Spec{K,L\S)) = Spec{K,H-'/\S)) = {{zi)}, (39) 

and for z ^ (zj) the unique solution pz G L'^{S) of {K — z)p = g is given by 

, 40 

. Zi- z 

with convergence in H~^''^{S). The corresponding formula for the scaled polarizability is 

«(-) = E;^' (41) 



where the sum is absolutely convergent. Finally, it holds that ^^UiVi = 2. 
Proof. Equation (j35p and Theorem 15.31 show that 

S^ec{K,H-^/^{S)) = Spec(^,L2(5)) = Spec(i^, l2(5)). 
Formula (P0]l also follows from ([35]) and spectral decomposition of A, since 



/L"} /O /L/'J Aj 



From a{z) = {pz, ^)l2(5) we now deduce (jH]). In terms of the measure p, this expresses the fact 
that p = "^iUiViSzi where J^. is the Dirac delta at Zj. This makes it clear that we have already 
proven that Yli '^i'^i = 2 as a part of Theorem 15.21 D 
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Remark 5.3. The statement that the -ff~"'^'^(S')-spectruin of K is equal to its L^(S')-spectrum is 
markedly untrue when S is not C^. In fact, when S is a square in two dimensions the H~^''^[S)- 
spectrum of K is real, while the -L^(S')-spectrum extends into the complex plane. 

We will now briefly discuss the limiting process which occurs when (Sk) is a sequence of sinooth 
surfaces approximating a non-smooth surface S. For brevity we assume that S is the cube defined 
by maxi<j<rf |rj| = 1 in d dimensions, r = (ri, . . . , r^), and that Sk is the superellipsoid defined by 
Si=i k«l^ = 1- We shall also skip some rather laborious details. See [50| for detailed approximation 
arguments when S" is a general Lipschitz surface. 

For any notation introduced so far, we denote by a subscript k that it corresponds to the 
surface Sk rather than S. Via appropriately defined homeomorphisms of Sk onto S we may, in a 
bicontinuous way, consider K^ and K* to be operators on the same space H^''^{S). We choose the 
implied isomorphism between if^'^(S') and H^''^[Sk) so that it extends to a unitary map of Lp'{S) 
onto L?'{Sk)- Similar conventions will apply throughout what follows. Using the boundedness 
results collected in Section [H and adapting the arguments in the proof of [50, Theorem 3.1], one 
can verify that K"^ converges strongly to K* on H^''^{S), meaning that liuikK^u = K*u in norm 
for any u G H^''^{S). Similarly, Kk converges strongly to K on H~^'^{S). 

Let z G C be a point at a distance at least e > away from the i/~^'^(S')-spectra of K and 
Kk, for all k. It is evident from pUj) that sup^ ||y0^^fc||_ff-i/2(5) < oo. We can hence extract a weakly 
convergent subsequence pz^k' with limit 6, 

Jii^(/5z,fe', w^>L2{s) = {b, w)l2(5), Mw G H^''^{S). 

Since s — \\ra.Kk = K we find that {Kk' — z)pz,k' converges weakly to {K — z)b. On the other hand, 
{Kk' — z)pz^k' = Qk' converges to g. Hence b = p^. Since every weakly convergent subsequence of 
Pz^k has the same limit pz, we conclude that Pz^k converges weakly to pz- In particular, noting that 
hk^hmH^/^{S), 

lim ak{z) = lim {pz,k,hk) l^is) = oi{z). 

Since this argument remains correct no matter the choice of gf G H~^''^{S) and h G H^''^{S), 
it is not difficult to see, based on the spectral representation (p6]) . that any z G Spec{K,H~^''^) 
is obtained as the limit of a sequence of eigenvalues Zk G Spec{Kk,H~^''^). See Weidmann [5^. 
From ([TT]) and ([5U|) we deduce that 

Spec{K,H'^/^{S)) = Spec{K*,H^/\S)) C [-1,1]. 

However, unlike the case of convergence in operator norm, we point out that strong convergence 
allows us to say very little about the character of the points in Spec(Er, H'^'"^). For example, they 
do not need to be eigenvalues. 
We summarize. 

Theorem 5.7. Let S he the unit cube in M , and let Sk be the approximating superellipsoid given 
by Ylii=i kil'^ = 1; ^ = (^ij • • • ^"^d)- Denoting by clR the closure o/i? C M, let 



B = cl 



[JSpec{Kk,H-'/''{Sk)) 
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Then Spec{K,H~^''^{S)) C B C [—1,1] and limfcQfc(z) = a{z) for all z ^ B. Furthermore, the 
supports of fik and /i are contained in B, and /ifc converges weak-star to jjl as measures on B. That 
is, 

lim / f{x) dMfe(x) = / f{x) d^l{x), V/ G C{B), (42) 

where C{B) denotes the continuous functions on B. 

Proof. Only the final point remains to be proven. However, the argument preceding the theorem 
can be repeated to show that (|42p holds for /(x) = \/{x — z)*, where z ^ B and t > is 
an integer. Now (|42p immediately follows in general from the fact that the linear span of such 
functions l/(x — z)* is dense in C{B), which in turn follows for example from the Stone- Weierstrass 
theorem. D 

We turn to the discussion of /i and a for a general Lipschitz surface S. First note that the 
support of // is contained in the iJ~^'^(5')-spectrum of K, and that fj, has a unique decomposition 

Here ^a is the absolutely continuous part of the measure, so that ^a = fJ-'ix) dx, where /u' G L^(E). 
fl has at most a countable number of atoms Zj, corresponding to bright plasmons at Zi, and fip is 
the atomic part of fi, 

i 

Note that each atom arises from an eigenvalue z of K, but not necessarily every eigenvalue is 
given positive measure by jjl, reflecting the distinction between bright and dark plasmons. Finally, 
Us denotes the singular part of //, excluding atoms. This means that /i^ has no atoms, yet lives 
solely on a set with Lebesgue measure zero. That is, there exists a measure zero set Rq such that 
fis{R) = for any Borel set i? C M with Rr\RQ = %. 

For < p < oo, let W^onf(C+) denote the conformally invariant Hardy space on the upper 
half-plane C-|- , consisting of analytic functions / in C+ such that 

ll/llp = sup/ f^TCdx<oo, z = x + iy. 

r>0 Jy=r F + *l 

This is known as the conformally invariant Hardy space as 

^conf(C+) = {/°a; : f G H^iB)}, 

where HP(0) is the usual Hardy space of the unit disk and a; is a conformal map of C+ onto the disk. 
'^conf(*^+) does not coincide with the usual Hardy space of the upper half-plane. Of importance 
to us is the fact that every / G ^conf (*^+) ^^^ boundary values /(x) = limy_^o+ fi^ + %) almost 
everywhere, and 



1 + x^ 



See Koosis [33] ^or further information. 

As previously mentioned in the proof of Theorem 15.21 a is the Cauchy transform of /i, a{z) = 
/C^(z). Changing variables in the Cauchy transform through u, exploiting that fi is compactly 
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supported and applying Smirnov's theorem, see [8|], leads to the fact that a G ^conf('^+) ^^'^ 
< p < 1. Hence, a has boundary values q^[x) for almost all x G M, taken as limits from the 
upper half-plane. A similar discussion could be carried out for the lower half-plane, but since ^i is 
real, the boundary values obtained from below are related to those obtained from above simply by 
conjugation. 

The absolutely continuous part of [i is related to a^ through equation ([25]) almost everywhere, 
7rfi'{x) = Q{a^{x)). Atoms Zi can be recognized as those points where \a{zi + ie)\ ~ £~^ as e — )• 0+. 
While some results are available, see [8,], recovering fis from a is much more subtle. However, if it 

turns out to be the case that 

|q:"'"(x)' 



/ 

Jr 



■dx<oo, (43) 

1 + x"' 

it follows that a G 7il^j^^{C+). But then / o u}~^ is the Cauchy integral of its boundary values and 
we infer that // must be absolutely continuous. That is, in this case ^p = //^ = and 



a[z) 



[iM±, z^S,eciK,H-^/\S)). (44) 

Jr X — z 



While we offer no strict proof, the numerical evidence in Sections [10] and [TT] suggests that ()43p 
holds when S" is a square in two dimensions or a cube in three, and hence that ()44p holds. 

To end this section we shall show that while a{z) has boundary values a'^{x) almost everywhere, 
the same cannot be said of the corresponding solutions pz- Note that the limit lim/^^^o /^ d/i 
exists finitely almost everywhere for x G M. We say that p'{x) exists whenever this is so, 

rx+h 

a'ix) = lim / dii, 

and as before we denote by a^ the set 

(7^ = {x G M : ^'(x) > exists}. 

Theorem 5.8. For x G M and e > 0, let Px+ie G H^^''^{S) denote the unique solution of 

{K - X- ie)px+ie = g- 

For any x G o"^, we have 

lim 1 1 Px+ie 1 1 H- 1/2 (5) = oo, (45) 

and moreover, the equation [K — x)p = g has no solution p G H^^'^{S). Hence, for such x, 
there exists no potential Ux solving the boundary value problem ^, ([2|), ([3]) in the sense given by 
Proposition I5.il 



Proof. If ([I5]) does not hold, there is a sequence (e^) with e„ — )• as n — >• c« and such that 
lim„_>.oo Px+i£„ exists weakly. From the proof of Theorem l5.2l we see that T^+i^^ = {K* — x — ien)~^h 
converges weakly in H^''^{S)/C to some element t^, meaning that 

Then, clearly, it must be that 

{K* - x)t, = h. (46) 
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Let ^T-a; be the positive measure defined by ^^(t) = {£{t)Tx,Tx)fji/2/gyQ, where S is the spectral 
measure of K* on H^''^{S)/C Recalhng that fi{t) = 2\V\{£{t)h,h) ^^1/2^^)/,^, we can reformulate 
(jMD as 

^l{t) = 2\v\{t-xf^lrAt)■ 

We infer 

This implies that limy_^o+ -^^(^ + iy) = for the Poisson integral 

PAz) = -It. To-, — 2' z = x + iyei:+. 

TT J^{t- xy + y^ 

On the other hand, it is well known that limy_^o+ P^li^ + w) = f^-'i^) whenever the right hand side 
exists finitely. This contradicts the hypothesis that /u'(x) 7^ 0. 

Suppose that there exists a p G H~^'^{S) such that {K — x)p = g. Note that 



{A-x-ie)~'^{A-x)= / ^^—d£A{t), 

jRt — X — IS 



where £a is the spectral measure of A on Lp'{S). Hence the operator norms || (A — x — ie) ^{A — x)\\ 
are uniformly bounded for e > 0. Since 

Px+ie = (K — X — ie)~^{K — x)p = VS {A — x — ie)~^{A — x)\fSp 
we have obtained a contradiction to (1451). D 



Remark 5.4. As to be expected, a similar argument shows that the conclusions of Theorem 15.8 
are valid also when x G M is such that /i({x}) > 0, that is, when there is a bright plasmon at x. 

Nevertheless, x >-^ px for x G M can always be given an interpretation as a H^^''^ {S)-'valued 
distribution, for example in the following weak sense. Let u : R — >• C be a Co^-function, let 
w G H^''^{S), and consider the functional 



Observe that 



p~^{u,w)= lim^ l^u{x){px+is,w)L2(^s)<ix- 



u{x){px+ie,w)L2(^S)^x= I u{x) I - — ^'"' dx = - / ICu{t-ie)dpg^w{t), 

Jr jR't — X — le jjR 

where Pg^w{t) = {£A{t)\Sg, \/S w)i2(^gy Due to the high regularity of n, passing to the limit in 
the Cauchy transform /C„(t — ie) poses no problems, and we obtain 



p+(n,u;) = - / lCu{t)dpg^yj{t). 

Jr 

Furthermore, p"*" solves the equation [K — x)/?"*" = 5 in the sense that 



p'^{u{x),K*w) - p'^{xu{x),w) = {g,w)L2(^s) /_u(x)dx. 
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Of course, we could equally well have introduced a solution p 



p {u,w) = lim / u{x){pa;-ie,w)L2rs)dx. 

We point out also that Kp = (/)+ + p^)/2 and Qp = —i{p'^ — /o^)/2 give us two real- valued 
distributions satisfying {K — x)'Rp = g and {K — x)Qp = 0. 

6. Capacitance 

The electrical capacitance C of an isolated charged conducting body V can be defined as 
the ratio of its total charge to its constant potential U{r). The problem of calculating C for V 
being a (unit) cube has "long been considered one of the major unsolved problems of electrostatic 
theory" [46| and attracted interest by researchers in computational electromagnetics for half a 
century. See [3|, |29|, |43| for recent contributions along with reviews of previous work and tables 
of historical progress. The highest relative precision for C so far, 10~^, was achieved in 2010 by 
parallelizing a Monte Carlo method and running it on a PC cluster [29]. See also Table [TJ 

The problem of determining C can be modeled as an integral equation much in the same way 
as the problem of determining a{z) and we omit details. If one solves 

(I + K + Q) p{r) = 1 , (47) 

where K is as in (jlOp and Q is the surface integral operator 



Qp= p{r) dar , 
Js 



(48) 



then the (normalized) capacitance can be evaluated as 



where 



''=4rf7M' ''=''• <"'" 



U{r)= G{r,r')p{r')dar>. (50) 

Js 



1. Strategies for computing a"'" (a;) 

The difficulties with computing q+(x) when S has edges and corners relate to issues of stability 
and resolution. We know from Theorem 15.81 that the electrostatic problem does not have a finite 
energy solution C/(r) and that ([6]) does not have a solution p(r) G H^^''^{S) for z € a^. Close 
to o"^, stability problems can be expected. A computational mesh needs to be extremely refined 
(locally) in order to resolve U{r) or p{r) and the solver must have the capability of dealing with 
strongly singular solutions. 

One strategy for alleviating these problems is to round edges and corners so that S becomes 
smooth. Then U{r) has finite energy and p{r) G L'^{S), except for at z = Zi, and a{z) assumes the 
form (|18p . Some distance above the real axis, corresponding to bigger losses, a{z) may resemble 
a'^{x) and can be evaluated using commercial software. This is essentially the approach in 3d.l5ll. 
IsJ], where finite element solvers are chosen. 
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Figure 1: Meshes on the boundary S of a square. Left: the coarse mesh and the corner vertices Sk- Middle: a mesh 
that is refined risuh = 3 times. Right: local meshes Ma., Mh and A^c centered around the corner vertex S3. 



It is, however, possible and advantageous to take the hmit z ^ a^ in a{z) nuniericahy while 
letting S retain its sharp shape. The rounding of edges and corners, while smoothing solutions, 
introduces new length-scales which is an unnecessary complication. In fact, there has been an 
intense activity in the area of constructing numerical al gori thms for solving integral equations on 
non-smooth curves in recent years [J, la, la, |7|, |23|, |2J, |26|, |27||. See [25|, Section 1.3] for an overview 
and a comparison of various approaches. Sharp corners and other boundary singularities can be 
treated extremely efficiently using fast direct and fully automatic solvers and by taking advantage 
of asymptotic self-similarity. An algorithm to this effect for a quantity analogous to a{z) for squares 
in a periodic setting is assembled and tested thoroughly in [25[. This paper continues along the 
lines of that work. 



8. Algorithm for the square 



This section is a summary of results from Refs. [23l . |25| | applied to the solution of (jlOp for V 
being a square. We construct two meshes on 5 - a coarse mesh and a fine mesh. The coarse 
mesh has 16 quadrature panels. The fine mesh is constructed from the coarse mesh by rigub times 
subdividing the panels closest to each corner vertex Sk in a direction towards the vertex. See Fig. (H 

8.1. Preconditioning and discretization 

Let S"^ denote a segment of S covering the four panels on the coarse mesh that lie closest to 
the corner vertex Sk - two panels on each side of Sk- The S^ are disjoint and their union is S. Let 
K{r,r') denote the kernel of K in (|10p . Split K{r,r') into two functions 



K{r, r') = K*{r, r) + K°{r, r) , 



(51) 



where K*(r,r') is zero except for when r and r' simultaneously lie on the same 5^. In this latter 
case K°(r,r') is zero. 

The kernel split ()5ip corresponds to an operator split K = K* + K° where K° is a compact 
operator. Discretization of ([TO]) , using a Nystrom method based on composite 16-point Gauss- 
Legendre quadrature and a coarse or a fine mesh on S, leads to an equation of the form 



(I + AK* + AK°)p = Ag, 



(52) 



where I, K*, and K° are square matrices and p and g are columns vectors. The matrix K* assumes 
a block-diagonal structure since the St are disjoint. This will be important in what follows. 
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The change of variables 

p{r) = {I + XK*)-^ ~p{r) (53) 

makes (|52]) read 



I + AK°(I + AK*)-^jp = Ag. (54) 

This right preconditioned equation corresponds to the discretization of a Predholm second kind 
equation with a composed compact operator. The solution p is the discretization of a function 
that is piecewise smooth and can be resolved by piecewise polynomials. 

From now on we let subscripts fin and coa indicate what type of mesh is used for the discretiza- 
tion. The collection of discretization points on a mesh is called a grid. The number ngub is assumed 
to be high enough so that pg^ resolves p{r) to the precision sought in our computations. 

8. 2. Compression 

The following decomposed low-rank approximation of the discretization of K° on the fine mesh 
holds to very high accuracy: 

Kg„ Rs PK°oj^P^ . (55) 

Here Kg^^ is a (256 + 96nsub) x (256 + 96nsub) matrix, K"^^^ is a 256 x 256 matrix, P is a prolongation 
matrix from the coarse grid to the fine grid and 

P^ = WfinPW,-i, (56) 



where W is a diagonal matrix containing the quadrature weights of the discretization, see [2c 
Section 5]. Superscript T denotes the transpose. 

The relation (|55p has powerful consequences for computational efficiency in the context of 
solving (|54p. As we soon shall see, it allows us to compress that equation and obtain the accuracy 
offered by the fine grid while working chiefly on the coarse grid. Only (I + AK*)~ needs the fine 
grid for resolution. We introduce the compressed quadrature-weighted inverse 

R = P^(Ifin + AK^J-ip. (57) 

With (j57p . the discretization of (jlOp assumes the final form 

(Icoa + AK°„^R) p^„^ = Agcoa , (58) 

where all matrices are 256 x 256. For later reference we introduce 

Pcoa = Rpcoa (59) 



as a weight-corrected density |27|, Section 5] 



8.3. Recursion 

The compressed inverse R has a block diagonal structure with four identical 64 x 64 blocks 
Rfc associated with the vertices s^- The construction of a block R/j from the definition (j57p is 
costly when n^uh is large. Fortunately, this construction can be greatly sped up via a recursion. 
In general situations this recursion uses grids on hierarchies of local meshes, see [23|, Section 6] 
and [24, Section 5]. For wedge-like corners, thanks to scale invariance of K{r,r'), only two local 
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meshes A^b ^-nd A^c are needed, see right image of Fig. [TJ The recursion for block Ry^ assumes the 
form of a simple fixed-point iteration 

R^fc = Pw-bc (lF{R^ii) J + lb + AK^fc) Pbc, i = l,...,nsub, (60) 

where Hik = Rfc for i = Uguh- The quadrature weighted and unweighted prolongation matrices 
Piybc and Pbc act from a grid on a local mesh Mc to a grid on a local mesh A^b- The superscript 
o in (1601) has a similar meaning as in ()5ip and the operator F{-} expands a matrix by zero-padding, 
see [25|, Section 6]. 

The derivation of ()60p relies on a low-rank approximation similar to ()55p 

K°p.PabK£Pkb, (61) 

where Ka is a discretization of if on a multiply refined local mesh A^a- See J26l . Section 7] for 
details. Conceptually, one could think of ([60]) as a process on a multiply refined local mesh, going 
outwards from the vertex, where step i inverts and compresses contributions to R^ involving the 
outermost panels on level i. 

The number ngub needed for resolution of R^ may grow without bounds as z approaches cr^sq 
(in infinite precision arithmetic). In order to accelerate the recursion we activate a combination 



of numerical homotopy and Newton's method when deemed worthwhile, see |25j, Section 6]. The 
Newton iterations are continued until the relative update in Rj^ is smaller than lOOemach or a 
maximum number of 20 iterations is reached, which roughly corresponds to a local mesh Ala that 
is refined rigub ~ 2^^^ « 10^ times. 

8.4- Solution, post-processing and interpretations 

Once the 256 x 256 linear system ()58p is solved for Pcoa' various quantities of interest can be 
computed. For example, the polarizability P^ becomes 

a{z) = h^oaWcoaRpcoa > (62) 

where h is the discretization of h{r). Results produced in this way are extremely accurate and fully 
confirm 24 of the entries for a{z) with |2;| > 1 in Table 1 of [39]. The remaining 5 entries differ in 



the last digit. Section 8 of |25|] gives error estimates for results produced in periodic settings. 

The original density pg^ can be reconstructed from ^^.^a by, in a sense, running the recursion ([60]) 
backwards (inwards on a multiply refined local mesh) . If this process is interrupted part- way, one 
is left with a mix of discrete values of the original density p (on outer panels) and quantities which 
can be easily converted into discrete values of a weight-corrected density p (on the innermost 
panels). Details of the process are given in [231, Section 7]. Here it suffices to observe that there 
exists a rectangular matrix Y, say, whose action on p^^^ produces entries of pg^. 

Let Q denote a restriction matrix from the fine grid to the coarse grid. Then 

Pcoa = QYpcoa (63) 

and 

K°oaR.Pcoa = K°„aR (QY)"^ p^oa • (64) 

We see that the blocks of the block-diagonal matrix R(QY)^ have an interpretation as multi- 
plicative weight corrections needed if K°Qa is to act accurately on p^oa- Other useful interpretations 
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Figure 2: Meshes on the surface S of a cube. Left: the coarse mesh. Middle: a mesh that is refined risub 
Right: local meshes A^a, A^b and A^c centered around a corner vertex. Compare Fig. [l] 



3 times. 



of the matrices introduced include: The columns of Y are discrete basis functions for p{r) on the 
fine grid; The columns of R are discrete basis functions for p{r) on the coarse grid multiplied with 
quadrature weights; the rectangular matrix Y (QY)~ maps p^oa to pgn. These observations will 
be used in what follows. 

The machinery of Sections 18. 1118.21 and 18. 31 is useful in several ways. The preconditioning aspect 
of ()58p reduces numerical error and improves the convergence of iterative solvers. The compression 
aspect of (|58p saves degrees of freedom and makes the algorithm fast and memory efficient. The 
recursion (|60p resolves the singular nature of p{r) close to corner vertices in an automated fashion 
and provides efficient basis functions and quadrature weights contained in the matrix R. No 
asymptotic analysis is required - we simply use Gauss-Legendre quadrature on the coarse mesh 
and on the local meshes A^b and M.c on which the recursion ([60]) takes place. Most, but not all, 
of these features can be retained as we step up into three dimensions. 



9. Algorithm for the cube 

Our algorithm for the cube mimics that of Section[8l A problem, however, arises in the split ()5ip . 
Unlike the square, the cube has both sharp edges and sharp corners and it is not possible to identify 
suitably disjoint surface elements 5^ that allow for an operator split K = K* + K° where K° is a 
compact operator. Thus, one cannot construct a block-diagonal matrix R which contains weighted 
basis functions for p(r) that simultaneously resolve the singularities stemming from the edges and 
from the corners. We shall circumvent this difficulty by focusing solely on the cube corners and 
construct the coarse and the fine mesh with square quadrature panels according to Fig. [51 that is, 
in complete analogy with Fig. [TJ As to compensate for the lack of refinement towards the edges, 
the discretization of K* and K° will incorporate ready-made one-dimensional basis functions and 
weights constructed for a square with the same A according to Section [ 



9.1. Preconditioning and discretization 

This section is a counterpart to Section [8. 11 Let the eight corner vertices of the cube be denoted 
Sfc and introduce surface element S^ that cover the 12 quadrature panels on the coarse mesh that 
lie closest to Sfc. Furthermore, let the smaller surface elements 5"^* be such that they cover the 
three panels on the coarse mesh that lie closest to Sk- The S*^ are disjoint, their union is S, and the 
kernel split ([5ip results in an operator split where the part of K° which accounts for interaction 
between points in 5* \ S"^ and S^* is compact. This restricted compactness is sufficient for our 
purposes since mesh refinement only takes place on the 5"^*, see the middle image of Fig. [2j 
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We discretize pU|) and make the change of variables that leads to (|54p . The discretization of 
K proceeds in three steps. First, we use the Nystrom method applying tensor products of rip- 
point Gauss-Legendre quadrature formulas on all quadrature panels to get an initial K. Then, for 
columns of K acting on discrete densities on panel pairs neighboring a single cube edge, but not a 
corner, we correct the Gauss-Legendre weights in the direction perpendicular to the edge. These 
corrections are realized by multiplying submatrices of K with blocks of Rgq (QsqYgq)" , where 
subscript sq indicates the square. See Section 18.41 This is our basic discretization. The resulting 
K acts accurately on p in situations describing the convolution of K(r, r') with p(r') both for r' 
away from edges and corners and for r' close to an edge but away from corners and r away from 
r'. 

Lastly, we change blocks of K describing interaction between panel pairs neighboring each other 
on opposite sides of an edge. Such interaction requires special attention due to the non-smooth 
nature of K(r, r') close to edges. We use interpolatory quadrature based on polynomial basis 
functions in the direction parallel to the edge and on the basis functions of Ygq in the direction 
perpendicular to the edge. This gives our final K. The total number of discretization points on 
the coarse mesh is n = 96np. The fine mesh has (96 -|- 72nsub)'^p points. 

The choice of the columns of Ygq as discrete basis functions for p{r) in the direction perpendic- 
ular to edges should be asymptotically correct, assuming that the singularities p{r) are dominated 
by two-dimensional effects away from the corner. Still, these basis functions are not optimal and 
they are responsible for the slower rate of convergence that we shall see when solving (j58p for the 
cube, compared to when solving ([^5|) for the square. 

9. 2. Compression, recursion and post-processing 

The compression of (j54p for the cube is analogous to that for the square in Section 18.21 The 
only difference, apart from that various matrices have different sizes, is that W, which contains 
the quadrature weights of the basic discretization and enters into the definition of prolongation 
matrix P^y of ([55]) . is no longer diagonal. The blocks of Rgq (QsqYgq)" , used as multiplicative 
weight corrections, are generally full matrices. 

The recursion for R of the cube, from now on denoted Rcu, follows Section E3] exactly. Some 
A allow for a rapid convergence in ()60p . Other A, corresponding to z close to parts of cr^cu, require 
that we resort to Newton's method and numerical homotopy. Note that recursions are carried out 
twice in the scheme: both for Rgq, needed for the discretization, and for Rcu- 

The polarizability a{z) of the cube can be computed from (j62p once the 96np x 96np system ()58p 
is solved. The matrix Wcoa in (|62p is now weight-corrected as described in the first paragraph of 
this section. We use the GMRES iterative solver for (j58p and make some use of symmetry in order 
to reduce memory requirements. 

10. From circle to square 

In a series of numerical experiments we now study the spectrum of K and the polarizability 
a{z) for a surface S that is gradually transformed from smooth to non-smooth. Such a study is of 
interest for several reasons. Due to the difficulties associated with solving electrostatic problems 
on non-smooth domains, see Section [3 it is common to round sharp boundary features prior 



to discretization [30|, l5l|, l5j] . Numerical effects, similar to those caused by rounding, could also 



result from insufficient resolution 4J]. Furthermore, no edge or corner in real world physics is 
infinitely sharp and the degree of edge smoothness can be critical in the design of, for example, 
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superellipse |r r+|r r=1 



superellipse |rj +\rj =1 




Figure 3: The eigenvalues Zi oi K (locations of plasmons) for the superellipse varies with k. Left: Zi with 2 < fc < 10^®. 
The eigenvalue at —1, a dark plasmon, is omitted. Right: zoom of positive eigenvalues with 2 < fc < 10. 



nanoantennas [19] . Finally, the experiments illustrate the theory overview of Section [3] in a setting 
which allows for high accuracy. All experiments are executed on a workstation equipped with an 
IntelCore2 Duo E8400 CPU at 3.00 GHz and 4 GB of memory. 
Similar to Klimov [32i] we let S be the superellipse 



I I h 

|n| + 



r2\ 



1 



(65) 



which for A; = 2 is a circle and for A; — )• oo approaches a square. We first compute eigenvalues of 
K using a discretization based on composite 16-point Gauss-Legendre quadrature and adaptive 
mesh refinement. Particular care is taken in the parameterization of S as to allow for resolution at 
boundary portions of high curvature. The accuracy in these computations varies with k. A rough 
estimate is a relative error of logxo(A;) • emach- 

Eigenvalues and the nature of their corresponding plasmons are shown in Fig. [3l The only 
bright plasmon at A; = 2 is a dipole. When k > 2, the bright plasmons have potential fields 
that are a mix of modes (dipoles, octupoles, etc.). The left image of Fig. [3] shows that K of the 
superellipse at A; = 10^^, which is close to a square in double precision arithmetic, has a spectrum 
that does not look continuous to the eye. The right image of Fig. [3] zooms in on the spectrum 
at low k. Klimov [32], in an analogous study for a superellipsoid with 2 < A; < 6, observes a 
phase-transition at A: = 2.5 and a critical point at A; ~ 3. 

The number of eigenvalues Zi with comparatively large residues UiVi in a{z) of (jlSp depends 
on A;. Fig. H] shows residues and polarizabilities a{z) for k = 10^^. Fig. [5] also shows how a{z) 
approaches a slowly varying function as z migrates from WJ^ = [—0.5, 0.5]. 

Fig. [5] compares a~^{x) with a{x + 0.05i) for a square. The algorithm of Section [8] is used. 
Each data point takes only a fraction of a second to compute and is accurate almost to machine 
precision except for x very close to {—0.5,0,0.5}. For example, the relative difference between the 



computed value of a^(l) and its known value — r(j) /Svr , see [49|], is 2 • 10~ . The left image 
of Fig. shows that the square has no bright plasmons (no poles in a~^{x)) and stands in forceful 
contrast the top right image of Fig. [U which exhibits a myriad of plasmons for the superellipse at 

k = 1012. 

It is interesting to compare the right image of Fig. [5] with the lower right image of Fig. HI Already 
at a distance of 0.05 away from the real axis, a{z) of the square and a{z) of the superellipse at 
k = lO^^^ are similar, as to be expected in view of Theorem 15.71 The convergence as A; — )• oo is 
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Figure 4: The number of eigenvalues Zi with comparatively large residues UiVi for the superellipse (dominant plasmon 
modes) depends on k, and a{z) of (I18|l decays away from z G [—0.5,0.5]. Here k — 10^^. Upper left: the 208 largest 
residues. Upper right: Qf(a;) with x £ [—1,1]. Lower left: Qf(a; + 0.01i). Lower right: q:(x -I- 0.05i). Compare Figs. [3] 
andO 



The square: polarizability 




The square: polarizability 




Figure 5: Polarizability of a square. Left: a"*" (a;). The curves are supported by 16492 data points whose relative 
accuracies range from machine precision to five digits. No convergent results were obtained within a distance of 10~^ 
from a; = 0. The values of ^{a+{x)} at a; = ^0.5 are ±10.3121215292. The sum rule (HU, evaluated via ^ using 
a composite trapezoidal rule, holds to a relative precision of 10^^. Right: a(x + 0.05i). 



uniform in z in compact sets away from [—1,1]. The accuracy achieved and the time required 
to evaluate a{z), however, are very different. While the accuracy in a{z) of the superelhpse, 
computed via (jlSp . is perhaps four digits and involves the eigenvalue decomposition of a 5376 x 5376 
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Figure 6: Convergence with the number of discretization points for the polarizability q+ (a;) and capacitance C of a 
unit cube. The values of x correspond to: cube with infinite permittivity {x — —1), resonance in corners (x = —0.6), 
resonance along edges (x — 0.25), and cube with zero permittivity {x = 1). See Table [l] for reference values used in 
error estimates. 



matrix, the accuracy in a{z) of the square, computed via ()58p . is much higher and involves only 
computations with matrices of size 256 x 256. We conclude that even when corners need not be 
strictly sharp from a physical viewpoint, it pays off to keep them sharp from a numerical viewpoint. 



11. The cube 

This section presents numerical results for the polarizability a^{x) and capacitance C of the 
unit cube produced by the algorithm of Section [9l While computing C has a long history in 
computational electromagnetics [3, l29|, |43|] , computing a+(x) is less explored territory [43] that only 
recently, with the rapid growth of the field of nanotechnology, has become fashionable. For example, 
•^/icu = {x ■ /^'(x) > 0} seems to be largely unknown. Fuchs [14]] and Langbein [30], over thirty years 
ago, found six or eleven approximate eigenvalues for K of the cube ("major/dipole absorption 
peaks") in the intervals [—0.573,0.408] and [—0.586,0.274], respectively and which supposedly 
account for around 95% of the sum ()19p . In view of the findings of the present paper, so far, one 
could suspect that these peaks are an artifact of insufficient resolution or unintended rounding of 
edges. Be as it may, this pioneering and highly cited work is now of interest in nanoplasmonics where 
the coupling of plasmon modes in metallic nanostructures such as nanocube dimers is important 
and interpretations seem to rely on these peaks [19|, l3^, l47l.l54]. 

In our algorithm for a~^(x) via ()58p . a complex valued R generally implies a complex valued 
solution Pcoa! complex valued basis functions Y, and 9{a+(x)} > 0. The square exhibits complex 
valued limits of Rsq whenever z = x + iy £ C+ -^ x £ cr^sq- The cube, which uses Rgq and Ygq 
for the discretization of K and which in turn enters into the recursion ()60p for Rcuj will therefore 



28 



have complex valued limits of Rcu and 9{a^(x)} > whenever x G (T^sq- We refer to this as 
resonance along edges. The reason being that 9{a^(x)} > implies x G Spec{K,H~^''^) and, as 
a consequence of the symmetrization arguments in Section [5] and the fact that the spectrum of a 
self-adjoint operator consists of eigenvalues and approximate eigenvalues, each such x is either an 
eigenvalue or an approximate eigenvalue of K on H~^''^{S). Moreover, Rcu may remain complex 
valued throughout the numerical homotopy process as z = x + iy £ C+ — )• x for other x S [—1,1] 
as well, that is, where limits of Rgq are real valued. We refer to this as resonance in corners. 

Fig. [6] shows results from convergence tests. The total number of discretization points on the 
cube surface is n. The examples with n > 2 ■ 10^ were executed on a workstation equipped with 
an IntelXeon E5430 CPU at 2.66 GHz and 32 GB of memory. Out of the chosen values of x, one 
can see that the error in a'^{x) is largest for x = 0.25. This is not surprising. A weakness in 
our algorithm is the assumption that the columns of Ygq are efficient basis functions for p{r) in 
the direction perpendicular to an edge. When x = 0.25, there is resonance along edges and the 
shortcomings of this assumption should be particularly visible. When x = — 1, x = 1 and for C 
there are real valued solutions p{r) S L'^{S) which are easier to resolve. 

Table 1: Reference values, estimated relative errors, and best previous estimates for the limit polarizability a'^{x) 
and the capacitance C of the cube. 
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previous results 


Ref. 
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a+(-l) 


3.644305190268 


10-^^ 


3.6442 


[48] 


3 • 10-^ 


a+(-0.6) 


5.85574775 + 16.642056431 


10-8 








a+(0.25) 


-2.76289925 + 3.080340351 


10-7 








a+(l) 


-1.638415712936517 


10-14 


-1.6383 


[48] 


6 • 10-^ 


C 


0.66067815409957 


10-13 


0.66067813 


[29] 


10-7 



The recursion (|60p for Rcu requires a substantial amount of memory when rip is large and 
Newton's method is activated. This explains why the test series for a~''(— 0.6) in Fig. [6] had to be 
interrupted at n = 69984, that is, for rip = 27. Timings vary greatly with x, rip, and the error 
tolerances that are set in recursions and iterative solvers. For safety, we set tolerances low and 
quote the following approximate computing times for n = 9600: a+(— 1), a"'"(l), and C took one 
minute each; a'''(0.25) took two minutes; a~''(— 0.6) took ten minutes. The test series for a'^(— 1), 
a^(l) and C all confirm previous available results, see TablelH and improve on these with between 
six and nine digits. 

Fig. [7] is our main numerical result. It shows that the isolated cube has no bright plasmons (no 
poles in a'^{x)). For validation, see the figure caption, we have used the two sum rules ([22]) . ([23]) 
and a third sum rule 



x'^n'{x) dx = 0.433464767896 , 



(66) 



which has been determined numerically for the cube 22] ■ Fig. [7] also shows that u^cu is approx- 
imately equal to the interval (—0.694526,0.5), possibly punctured at {—0.5,0}, and raises the 
intriguing question of what surface shapes S lead to connected a^,. Our computed (T„cu is broader 



than those estimated by earlier investigators [IJ, |32, l36[ using other techniques. We conclude, 



again, that working with sharp edges and corners is an advantage. 
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The cube: polarizability 




real{a"^(x)} 
imag{a"'(x)} 



Figure 7: Polarizability of a cube. The curves are supported by 1195 data points wiiose relative accuracies range 
from ten to five digits. No convergent results were obtained for q"'"(x) within a distance of 10"^^° from x = —0.5 and 
within a distance of 10~^ from x = 0. The maximum value of K{a^(x)} is 20.00826 and occurs at x ~ —0.694526. 
The sum rules (|22|l . (|23p . and (|66p . evaluated via (|25p using a composite trapezoidal rule, hold to a relative precision 
of 3 ■ 10~^ 4 ■ 10"^ and 6 • 10"^ respectively. 

12. Conclusions and outlook 

We have constructed an electrostatic solver for cube-shaped domains and used it to produce 
new results for some canonical problems. A particular characteristic of the solver is that it takes 
advantage of sharp edges and corners, rather than being a victim of them. 

A mathematical framework capturing the symmetric features of the double layer potential and 
its adjoint has been identified, allowing for analysis of the polarizability a{z) through spectral 
theory for self-adjoint operators. Moreover, this framework and its corresponding machinery is 
natural and satisfactory also from a physical viewpoint, as the potentials produced to solve the 
electrostatic problem are exactly those of finite energy. With full mathematical rigor we have shown 
how to, with respect to polarizabilities, interpret the limit process which occurs when deforming a 
smooth surface into a cube, a point which has generated a fair amount of discussion in the materials 
science community. 

Furthermore, we have shown that while the polarizability via limits can be extended to a 
function a'^{x) defined almost everywhere on the real axis, the same statement fails for the cor- 
responding potentials whenever the representing measure ^ has a non-zero absolutely continuous 
part. That is, no finite energy potential can be associated with a point j; £ M where n'{x) > exists 
and is non-zero, neither as a direct solution to the electrostatic problem, nor as a limit from the 
upper half-plane. We remark, however, that while the density distributions px hence fail to exist 
for such X, the whole map x ^ px can still be interpreted naturally as a vector-valued distribution. 

For the cube, the mathematical theory in conjunction with the numerical findings show that 
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d/i(x) = fj,'{x) dx is purely absolutely continuous, and that the set where fi'{x) > is given by the 
interval (a, b) ~ (—0.694526, 0.5), possibly excepting the two points x = —0.5 and x = 0. Hence the 
integral representation (j44p for a{z), in terms of /^'(x), holds. Finally, ^'(x) has been determined 
numerically. These discoveries, we hope, will be of use in plasmonics where absorption peaks of 
nanocubes play an important role. 

Future efforts will be directed towards solving the Helmholtz equation, compare [J,l5|,|6[. Should 
the polarizability of clusters of cubes or the effective permittivity of cubes in periodic arrangements 
be of interest, our solver requires only minor modifications. 
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